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Abstract 


Ground tests of solid propellant rocket motors have shown that metal-containing 
propellants produce various amounts of slag (primarily aluminum oxide) which is 
trapped in the motor case, causing a loss of specific impulse. Although not yet 
definitely established, the presence of a liquid pool of slag also may contribute to 
nutationai instabilities that have been observed with certain spin-stabilized, upper- 
stage vehicles. Because of the rocket’s axial acceleration-absent in the ground tests— 
estimates of in-flight slag mass have been very uncertain. Yet such estimates are 
needed to determine the magnitude of the control authority of the systems required 
for eliminating the instability. A test rig with an eccentrically mounted 
hemispherical bowl was designed and built which incorporates a "follower" force 
that properly aligns the thrust vector along the axis of spin. A program that 
computes the motion of a point mass in the spinning and precessing bowl was 
written. Using various RPMs, friction factors, and initial starting conditions, plots 
were generated showing the trace of the point mass around the inside of the fuel 
tank. The apparatus will be used extensively during the 1990-1991 academic year 
and incorporate future design features such as a variable nutation angle and a film 
height measuring instrument. Data obtained on the nutationai instability 
characteristics will be used to determine order of magnitude estimates of control 
authority needed to minimize the sloshing effect. 


Introduction 

Many rocket motor solid propellants in current use contain a significant amount of 
aluminum which, when burned, produces a slag consisting of aluminum oxide and 
elemental aluminum. Most of this material is expelled throughout the rocket 
motor nozzle and adds to the thrust, but some remains trapped in the motor case. 
The melting point of the (X-form of AI 2 O 3 is about 2050°C, below the temperature of 
the combustion gas. The liquid slag, in the form of small droplets, is subject to a 
combination of forces that include the drag from the combustion gas, the inertial 
force resulting from the axial acceleration of the rocket, and (for spin-stabilized 
vehicles) the centrifugal force resulting from the vehicle spin. 

The present analysis postulates that, because of the high level of turbulence in the 
motor, slag droplets entering the gas stream are ejected, and that trapped slag is 
formed primarily by liquid slag flowing along the surfaces toward the point of 
minimum potential energy in the accelerating and spinning motor. Also, the 
present analysis concludes that slag will accumulate to some degree in all spinning 
or accelerating rocket motors with aluminum-contaning propellants and submerged 
nozzles. 
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A number of spin-stabilized vehicles that use aluminized propellant have shown a 
marked tendency for a "coning” instability; i.e., a precession with steadily increasing 
nutation angle. These motors have a submerged nozzle geometry, resulting in a 
downstream annular pocket which is likely to favor slag retention. It has been 
surmised, therefore, that the sloshing motion of a liquid slag pool may be a 
contributing cause of the observed flight instability. The effects of liquid slag on the 
stability of spinning vehicles is similar to the effects produced by fuel slosh in 
spacecraft. Slag retention also requires examination because of its potentially 
deleterious effect on specific impulse. 

Through installation of witness plates downstream of the nozzle, where some of the 
(now solid) slag particles are deposited, estimates of the size distribution and total 
mass of the expelled particles have been made. Ground tests of this type, however, 
take no account of the processing of the droplets in the nozzle. 

This report consists of a mechanical design that simulates the motion of a spherical 
fuel tank in a thrusting spacecraft. A true simulation of the thrust was thought to be 
impossible due to the gravitational support forces present in the laboratory. 
However, through the means of an eccentrically mounted spacecraft model on the 
top of a turntable, the simulation of thrust aligned with the vehicle axis is possible. 
The mechanical design was finished during the 1990 winter quarter and the test rig 
was built in the spring. The comparison of the initial description (see figure 1) with 
the design actually built (see figure 2) shows the evolution of the design concept. 
Qualitative analysis will be provided by photographs of fluid profiles at given time 
intervals and quantitative analysis by correlation of film thickness from capacitance 
measurements between two platinum wires located in the bowl. This sensor will be 
designed, built, and incorporated into the test rig slip-ring assembly during the 1990- 
1990 academic year. From these data, nutational instability characteristics and order 
of magnitude estimates of control authority needed to eliminate the instability will 
be determined. 

A computer program was written to simulate the shape of a fluid in a spinning and 
precessing container with a nutation angle equal to zero. The fluid was assumed to 
be in hydrostatic equilibrium. The fluid depth as a function of position along with 
the shoreline of the fluid was determined. A more general code was written which 
computes the motion of a point mass in a spinning and precessing hemispherical 
container. Using various RPMs and friction factors, plots were generated to 
compare the motion of the point mass and validate the theoretical model (see 

figure 4). 
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1. general bearing assembly 

2. AC motor (variable rpm) 

3. pulley for motor shaft 

4. main drive pulley 

5. secondary pulley (stationary) 

6. main shaft 

7. control arm 


8. bowl pulley 

9. bowl bearing housing assembly 

10. bowl mounting flange 

11. hemispherical bowl (lucite) 

12. bowl suport shaft 

13. idler guide 


Figure 1: Appratus Diagram (Not to scale) 
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Figure 2: Completed Test Rig 

(a) top view showing liquid sloshing in bowl 

(b) front view showing dual motor assembly 
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Viscous Dissipation: 

The degree of instability of a thrusting, spin stabilized spacecraft depends strongly on 
the amount of internal energy dissipation. The dominant energy dissipation 
mechanism is thought to be caused by the sloshing of liquid slag at the bottom of the 
solid motor casing which directly influences the body’s motion. Oscillatory, and 
sometimes violent, motion of the fluid induce corresponding oscillations of the 
body Viscous effects in the fluid also influence the body causing the nutation angle 
to change thereby affecting the stability. It is, therefore, important to estimate the 
energy losses in the fluid. 

Once these energy losses are estimated, one can predict the body motion by reducing 
its kinetic energy at the same rate. This approach is known as the "energy sink" 
procedure. Due to the growing nutation angle from energy dissipation, thrust 
corrections need to be fired to stabilize the craft. This requires more fuel to be 
included for stabilization during launch which ultimately increases launch mass. 
Having to fire these correcting thrusters at the right time creates yet another 
problem in the attitude dynamics and control of the spacecraft. Ideally, nutational 
instability characteristics and order of magnitude estimates of control authority 
needed to eliminate the instability would allow designers to provide the lightest 
control system necessary to minimize this phenomenon. 


Scale Model Principles 

Many different models have been developed to test sloshing and its effect on 
spacecraft. Most of these models, however, are made to simulate the sloshing of a 
spacecraft in which thrust is absent. One of the recent problems is that an instability 
evidenced by a growing nutation angle has been observed during the firing of liquid 
and solid perigee and apogee motors. A new model to simulate this motion was 
needed which properly aligns the "thrust" vector with the model axis. 

A simple design of a spacecraft model mounted eccentrically on a turntable can be 
used. This rig simulates the thrust as a "follower" force (see figure 3). Previous 
models were subjected to gravity forces acting at the center of mass. But the new 
model produces a resultant of gravity and inertial forces that remains aligned at all 
times with the vehicle axis. Hence, this thrust "follows" the model as is spins and 
precesses around on the turntable. 
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Figure 3: "Follow Force Diagram" 


10 Second Marble Trace 

40 RPM, friction factor = 0.5 
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Because of space and cost constraints, it is necessary to have a model that is not full 
scale. It must then be shown that the model behaves in the same way as the 
spacecraft. Therefore, it is required for the model to have the same inertia ratio as 
the spacecraft: 


riii = m 

L^PJmodel L^PJspacecraft 

It follows that the ratio of the precession rate to the spin rate be the same in 
both the model and the full scale model. To simulate the dynamics of the sloshing 
requires that the Froude numbers of the model and spacecraft be the same: 


Froude number 


RtfdO/dt) 2 " 

RKdfc/dt) 2 ] 

. g/ COS0 _ model 

L T/M ! 


Solving for (dO/dt) model : 

MO/dO^e. - 

Using these equations, a good approximation to a thrusting spacecraft can be made 
in the laboratory. 


Mechanical Design: 

A distinct design evolution was experienced in attempting to construct a test rig 
which would adequately simulate the conditions present during the burn of a solid 
propellant rocket motor. As a preliminary experiment it was primarily designed to 
provide a qualitative analysis of fuel and slag sloshing and aid in the development 
of future experimentation. 

The design problem was to simulate rotation about the rocket’s own axis and the 
subsequent precession about an associated axis, both of which are effects of spin 
stabilization. It was initially agreed that dual rotating shafts were best fitted to 
produce the kinematics of the situation, and subsequently the design problem was 
limited to developing a system that would drive the two shafts with correct 
direction and rates of spin. In order to achieve this effect several proposals were 
made, first of which entailed using a set of belts and pulleys driven by a single 
electric motor. Succeeding designs included such elements as a planetary gear 
system, a set of rubber wheels, or a set of dual motors. In the end, the initial concept 
of belts and pulleys was adopted for their availability and ease of use. 
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The rig is mounted on a half inch thick aluminum table, approximately one meter 
square and held up by four nine inch long aluminum legs. The main shaft is 
positioned vertically through the middle of the table, housed by a bearing assembly 
which is mounted to the under face of the table. This shaft is driven by a belt, 
connected to a variable RPM electric motor also mounted from beneath the table. 
To the top of the main shaft is mounted a control arm made from an aluminum T 
beam. On one side of the control arm is the fuel tank assembly and on the other, an 
equal counter weight made of lead plates. 

The hemispherical bowl, turned from a Ludte block, is mounted to a second shaft 
which rotates within the bearing housing mounted to the control arm. Positioned 
on the main shaft and on the bottom of the second shaft are two pulleys. The pulley 
on the main shaft is secured and remains stationary with respect to the table. The 
other pulley is secured to the second shaft and produces the rotation of the bowl 
about its own axis. A crossing belt connects the two pulleys, and as the main shaft 
rotates at an average rate of forty RPM, the second shaft rotates twice as fast in the 
opposite direction. In order to keep an adequate tension in the belt, the bearing 
assembly housing the shaft, can shift horizontally by ± 0.5 inch. In addition, an idler 
is included on the control arm to guide the belt and maintain its tension. 

During the next academic year (1990-1991) a sensor will be designed which 
determines the film thickness by measuring the capacitance between two platinum 
wires. This will hopefully provide a means to quantify the force and momentum 
produced by the rotating liquid in the bowl at various RPMs. In order to incorporate 
this instrument, an electric connection to the bowl is needed through a set of slip 
rings in the rotating mechanism. Just below the bowl and above the bearing 
assembly is mounted the first slip ring. And at the bottom of the main shaft below 
the bearing assembly is mounted the second slip ring. To connect the wires from 
the control arm to the second ring, a hole is drilled down the center and through the 
entire length of the shaft. Through this hole the wires are run to the slip ring. 


Computer Simulation: 

A theoretical analysis which approximates the fluid in the bowl with a point mass 
was developed. The result was a system of two ordinary differential equations 
which can be solved numerically by Heim's method for initial value problems. A 
code was generated which determines the x, y, and z coordinates of a "marble" 
rolling around inside the bowl given a friction factor, initial starting coordinates, 
bowl RPM, and nutation angle. The friction factor was varied to simulate the effects 
of fluid viscosity and friction of the point mass. The larger the friction value, the 
more of a damping effect the marble exhibited. For smaller values, the marble took 
longer to stabilize and rose higher in the bowl (see figures). When the actual 
experiments begin this fall, the code can be properly validated with better estimates 
of the friction factor, RPM, and nutation angles necessary to demonstrate a valid 
theoretical model and test rig. 


9 


7t 



ORIGINAL PAGE IS 
OF POOR QUALITY 


TABLE OF CONTENTS 


INTRODUCTION 

LIST OF SYMBOLS 

REFERENCE FRAMES ■ • 

STABILIZATION PRINCIPLES 

SCALE MODEL PRINCIPLES 

MECHANICAL DESIGN 

DESIGN SKETCHES 

COMPUTER ANALYSES - A. SIMULATION OF -LUID SHAPE 

B. SIMULATION OF POINT MASS MOTION 

HOT TUNGSTEN WIRE ANALYSIS 

LIST OF MATERIALS 

CONCLUSION 

APPENDIX A: PRELIMINARY DESIGN 

APPENDIX 3: CALCULATION FO FLUID SHAPE SIMULATION 

APPENDIX Cs POINT MASS MOTION HANDOUT 

H ppcrf,jDix D: POINT ‘IASS MOTION COMPUTER PROGRAM LISTING 

APPENDIX E: VELOCITY CALCULATIONS FOR HOT WIRE 

REFERENCES 


I 

3 > 

H 

S' 

s 


13 

16 


HZ 

HI 

M\ 

.6fT 

Db 

.?o 

.sr 


10C 

106 


BIBLIOGRAPHY 


ORIGINAL PAGE IS 
OF POOR QUALITY 


INTRODUCTION 

In theory, the definition of a. rigid body does not 
permit any energy dissipation. However, it is known that 
all spacecraft and rockets have some non-rigid entities 
including elastic structural deflection and the liquid 
motion of fuel in its tanks, otherwise known as slosh. 
Since Explorer I, it has been known that these properties 
can have a major effect on the motion, that is, there can be 
instabilities that depend largely on the internal energy 
dissipation. In most spinning spacecraft and rockets, the 
largest amount of energy dissipation comes from the liquid 
si osh . 

This report consists of a me c ran design that 

simulates the motion of a spherical fuel tank in a thrusting 
spacecraft. A true simulation of the thrust was thought to 
oe impossible due to the gravitaicnal support forces present 
n the ’aboratory. However, thrcugn the means of an 
eccen tr i cal I y mounted spacecraft model on the too of a 
turntable, it was discovered in reference 1, *hat the 
simulation of thrust aligned with the vehicle axis was 
ndeeo possible. The mecnan i cal design is scheduled to be 
built and tested during the ~ ? " - e su . 1 Ue 

in the ^orm of pictures and data collected about the depth 
of water at pertain points. The results will then be 
pua 1 i tat i ve 1 y presented with conclusions drawn about the 
iiosnirg or the -usi . H cacKground on tr.e design ana a : ist 
pf mater a Is needec for the -design are ncluded. 

\ 


Also presented in this report is a theory on the 
possibility of using a heated tungsten wire to measure the 
depth of 1 iquid in the model . By relating the heat loss of 
the wire to the depth of fluid inside the cup, it was 
thought that an accurate reading of the depth could be 
obtained. However, after many calculations it was found 
that there were too many uncertainties for the hot wire to 
be an accurate device. Nevertheless, these calculations are 
valuable and therefore included. 

A computer program was written to simulate the shape of 
a fluid in a spinning and precessing container with a 
nutation angle equal to zero. The fluid was assumed to be 
in hydrostatic equilibri urn . " r he fluid depth as a function 
of position along with the shoreline of the fluid were 
determined. A program that determined the motion of a point 
mass in a spinning and precessing hemispherical container 
was also written. Using different conditions, • . e? . , 
different RPMs and friction factors, plots were generated to 
compare the motions of the point mass. 
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LIST OF SYMBOLS 


<X,Y,Z> 

<x,y,z) 

♦ 

e 

M 

m 

T) 


r 

a 

P 

9 

I 

T 


T9 



h 


D r 


Nu 


Inertial reference frame 

Body-fixed frame 

Precession rate 

Spin rate 

Nutation angl e 

Total mass including fuel 

Fluid mass 

Ki nemat i c v i sc os i ty 

Rad i u s 

Centripetal acce 1 er at i on 
Dens i t y 

Acceleration of gravity 
Momen t of inertia 
Thr ust 

Reynol ds number 
-roude number 
Pertaining to spacecraft 
Pertaining to model 
Height 
D r essur e 

Angular velocity t subscripts cbv i cus) 
Prandt! number 
Sr asftof number 
Nusse 1 t number 
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STABILIZATION PRINCIPLES 


There are primarily two different stabilization 
techniques currently in use. They are spin stabilization 
and three-axis stabilization, each of which has its 
advantages and disadvantages. Since we are concerned with 
rockets and thrusting spacecraft, we will be most interested 
in spin stabilization. 

Spin stabilization is based on the gyroscopic stiffness 
produced due to the rotation of the rocket or spacecraft. 
Dual-spin, which involves two types of spinning and also 
precessing of the rocket is the most common stabilizaion 
technique as opposed to single-spin. Since spacecraft and 
rockets carry much of their fuel in the form of a liquid, 
the motion of this fuel due to the spin stabilization and 
the presence of thrust is of great interest. This motion 
- r om its n cm i n a 1 position ' s wh at is K n own as 11 1 i q u i d 
s 1 osh i ng“ . 

— 0 o s c» r e e c +• • -• « ^ ^ ^ ; tv of t r. ? = o i nr < r - c, i-ccy c e d e n d s 
strong’ y on the amour.* :f internal energy dissipation. The 
most dominant energy dissipation mechanism s the -fuel 
si osh i no which ; nfluences the body's motion immediately and 
directly. Oscillatory, ana sometimes violent, motion of the 


fluid 

nduce 

correspond! ng esc 

illations 

of the 

body , 

V> i scous 

effect 

s in the f 1 u i c al so 

; n f 1 u e n c e 

the bod y 

c au sing 

the nut 

a t I on 

angle to change 

t h e r e by 

affect! 

n g the 

stab i 1 i t 

y . I 

t is j therefore 

mpor tan t 

to e s t i m 

ate the 

energy 

- 550 = 

~ c e - ■ w i c . 
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The energy dissipation *"■*. te :s assumed to be a 
of the many spacecraft parameters as follows, 

where , 


M = total mass including fuel 

m = fluid mass 

l) - kinematic viscosity 

R ■ radial distance of tank center from spin axis 
= radius of tank 
= radius to free surface 
a = centripetal acceleration, r < * > 

etc. 

expanding in an infinite series representation, 


<90 


t'l (K; d'-fll 

I s I 


M (f «' %L ) 


, < 4 o ; — t "? OO 

iAl h ere e i i = - x P on e n * s * - - 1 , 2 , . . . l * « ■ - 1 • 


the 


Usino the non-dimensional i zat : on process using 
1 Buckingham Pi ^heorem" (see reference 2) snd gr i t : ng the 
r p su 1 t in the original form (where terms of like exponents 
have been brought together and the exponents dropped), 



ome of these terms can be inverted to show 


their more 


- am i 1 i ar 



= Reynolds number 
= Fr oude number 
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Once these energy 


’ osses are estimated, one can then 


predict the body motion by reducing its kinetic energy at 
the same rate. This approach known as the "energy sink" 
procedure is beyond the scope of the class. A more detailed 
discussion of it and energy dissipation can be found in 
reference 3 and reference 4. 

One of the results of this growing nutation angle from 
energy dissipation is that correcting thrusts will have to 
be f i r ed to restab ilize the craft. Th i s w i 1 1 use up mor e 
fuel and shorten the 1 ife of the craft. Having to fire 
these correcting thrusts at the right time creates yet 
another problem n the a.tt: tube zyr am c= ^nc — n .. - • ■ 
The simplicity of spin stabilization means that ’iquid 
sloshing will still remain a challenging problem. 
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SCALE MODEL PRINCIPLES 


Man y ditterent 

mode 1 

s have been 

deve loped to 

test 

s 1 osh i hq and its 

effect 

on spacecr at t . 

Most of 

these 

models, however, are 

made 

to s i mu 1 ate 

the 

si osh i ng 

of a 


spacecr at t 

i n 

which 

thrust 

i s 

absent. One 

ot 

the recent 

probl ems i 

s that 

an 

i nstab i 1 

i ty 

e v 

i denced 

by 

a growing 

nutation 

angle 

has 

been observed 

during 

the 

firing of 


liquid perigee and apogee rocket motors reference 1). 
Therefore, a new model to simulate the motion had to be 
■found . 


As seen in reference 1, an intricate model need not be 
made. Instead, a rather simple design of a spacecraft model 
mounted eccentrically on a turntable can be used. This rig 
simulates the thrust as a "follower force" (figure 1). As 
opposed to figure 2, where the model is subject to gravity 
reacting forces acting at its center of mass, figure 3 shows 
the -esiqn jhere "he model lies eccentrically on a 
turntable. This produces a resultant of gravity and 
inertial -orces that remains aligned at all times with the 
vehicle axis. Hence, this thrust is of the type ■-ef erred to 
a "follower ^orcs" . 

Because of space s. no cost constraints, ,t ■ -- :ecessar> 


to have a model that is not full scale. It must then be 
shown that the model behaves : n the same way as the 
spacecraft. Therefore, it is required *or the model to have 
the same inertia ratio as the spacecraft, 
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It also -follows that the ratio of the precession rate to the 
spin rate be the same in both the model and the -full scale 
model. To simulate the dynamics of the sloshing requires 
that the Froude numbers of the model and spacecraft be the 


same 


F« 




(*T/w) 

From this it is seen that 


Uu 



thrust i ng 


o an be made . 
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MECHANICAL DESIGN 


A distinct design evolution was experienced in attempting to 
construct a test rig^ which would adequately simulate the 
conditions present during the burn of a liquid propellant rocket 
engine. As a preliminary experiment it was primarily designed to 
provide a qualitative analysis of fuel sloshing and aid in the 
development of future experimentation. 

The design problem was to simulate rotation about the rockets 
own axis and the subsequent precession about an associated axis, 
both of which are effects of spin stabilation. It was initially 
agreed that dual rotating shafts were best fitted to produce the 
kinematics of the situation, and subsequently the design problem 
was limitai to developing a system that would drive the two 
shafts with correct direction and rates of spin. In order to 
achieve this effect several proposals were made, first of which 
entailed using a set of belts and pulleys driven by a single 
electric motor. Succeeding designs included such elements as a 
planetary gear system, a set of rubber wheels, or a set of dual 
motors. In the end, the initial concept of belts and pulleys was 
adopted for their availability and ease of use. The details of 
test rig's design are the following: 

The rig is mounted on a half inch thick aluminum table, 
approximately one meter square and held up by four nine inch long 
aluminum legs. The main shaft is positioned vertically through 
the middle of the table, housed by a bearing assembly which is 
mounted to the under face of the table. This shaft is driven by 
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a belt, connected to a variable RPM electric motor also mounted 
from beneath the table. To the top of the main shaft is mounted a 
control arm made from an aluminum T beam. On one side of the 
control arm is the fuel tank assembly and on the other, an equal 
counter weight made of lead plates. 

The hemispherical fuel tank, turned from a lucite block, is 
mounted to a second shaft which rotates within the bearing 
housing mounted to the control arm. Positioned on the main shaft 
and on the bottom of the second shaft are two pulleys. The pulley 
on the main shaft is secured and remains stationary with respect 
to the table. The other pulley is secured to the second shaft and 
produces the rotation of the fuel tank about its own axis. A 
crossing belt connects the two pulleys, and as the main shaft 
rotates at an average rate of forty RPM, the second shaft rotates 
twice as fast in the opposite direction. In order to keep an 
adequate tension in the belt, the bearing assembly housing the 
shaft, can shift horizontal^ by + or - half an inch. In addition, 
an idler is included on the contol arm to guide the belt and 
maintain its tension. 

If a feasible means of determinir^ the depth of the fluid could 
be produced, an electric connection to the fuel tank would be 
necessary. For this purpose, incorporated into the design is a 
set of slip rings to make possible electric connections to the 
rotating mechanisms. Just below the fuel tank and above the 
bearing assembly is mounted the first slip ring. And at the 
bottom of the main shaft below the bearing assembly is mounted 

To connect the wires from the control arm 


the second slip ring. 


to the second ring, a hole is drilled down the center and through 
the entire length of the shaft. Through this hole the wires are 
run to the slip ring. 

The following are detailed drawings of the individual parts. 
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SCALE 1-1 



Parts 10&11. Hemispherical tank and flange 
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Simulation of Fluid Shape 


In order to gain insight into the instability of a 
spacecraft due to liquid slag, a simple study of a liquid 
under simulated static forces should be made. This study can 
be made by applying Bernoulli's Equation (1) for both 
spherical and Cartesian coordinate systems in 3 dimensions. 

-(V*P + ^"g^Vh) = ^ 1 ^ 

If applied properly a shoreline- -the line where the 
fluid meets the cup- -can be graphed to obtain a top view of 
the cup. Moreover a cross-section of the front view may also 
be obtained. Due to the difficulty that the actual spacecraft 
precession poses with its nutation angle, a simplified study 
shall be made using a cup without a nutation angle. 

Before the detailed study of the fluid shape under 
static forces is made, it is useful to conduct a qualitative 
examination of the fluid shape. Since the cup is spinning 
about its own axis, as well as about a point located at a 
fixed distance from the cup's center, one can assume that the 
fluid's actual shape is the superposition of the shapes 
caused by each motion separately. Figure 1 shows a picture of 
the posed problem while figures 2,3 and k are purely 
qualitative drawings of the fluids shape before any 
calculations were done. 


sr 







has been made, a theoretical solution should be constructed 
using equation (1) applied to our posed problem of fig. 1. 

Applying Bernoulli's equation in the x direction, eqn. 
(1) reduces to eqn. (2): 



( 2 ) 


where a„ is : 

3 m = Wrrup ^ *X ^ W m r m ^ * ( X* r m + X) Where • Ter up i X i Tc= up 


I > acceleration due to the cup's 

rotation about a fixed point. 

> acceleration due to rotation 

about the cup ' s axis . 

After integrating and assuming zero gauge pressure 
outside of the cup, equation (3) is obtained in the xy- 
plane: 


z(x) 


Wc lip ^ ^ Wmr- 


2 * g 


2 

-*X 2 


+ 


W* r IT, 2 • r- rrt 

* x + z ( 0 ) 


g 


(3) 


Similarly applying eqn. (1) in the y direction, the 
following is obtained (4): 



where a v is : 


( 4 ) 


ORIGINAL PAGE IS 
OF POOR QUALITY 



ay = Wc= UP . 2, y 


where : 


Tcup 


< y < rcup 


Note that there is only acceleration due to the cup's 
rotation about its own axis in the yz -plane at the instant 
time we apply eqn. (1). Following the same integrating 
procedure used to obtain equation ( 3 ) , equation ( 5 ) can 
easily be shown: 


Wr Lip 2 

z ( y ) = *y 2 + Z(0) 


Combining equations (3) and (5) we get a general three- 
dimensional expression for the height of the fluid in 
Cartesian and spherical coordinates (6a) & (6b). See appendix 
for the definition of spherical coordinates used. 


Wcup ^ + w*t-m 2 W«Lr-m 2 # r«n- 


z(x,y) = 


>x 2 + 


•X + 


2«g 


Wcr u p 

+ •y 2 + z(0,0) 

2 • g 


(6a: 


z(R,<t>,0) = 


W|-L.;ri 2 ^ Wj.r-,1 


2 • g 


-*R» cos ( <J> ) cos ( 0 ) + 


«R» cos ( <t> ) cos ( 0 ) + z(0,0) 


( 6b) 


In the xz -plane (front view, z(x,y=0) ), we obtain a 
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parabolic equation. By adjusting the constant z(0,0) one can 
determine how high the fluid may be at the middle of the cup 
before it would spill out. A quantitative plot was made and 
the initial height of the fluid determined to be 0.91cm using 
MathSofts's MathCad program. The following graph, constructed 
for the xz-plane, indicates the qualitative drawing for the 
final fluid shape in the xz-plane was indeed correct . 



To obtain a shoreline of the fluid ( xy-plane , top view) 
one must combine the equations for the height of the cup 
with the equation obtained for the height of the fluid shape 
(6b) . 

The since the cup is a sphere, the equation of the cup 
may be described by equation (7a) or (7b). 
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z 2 + x 2 + y 2 = R 2 


(7a) 


z = *J R 2 - x 2 - y 2 (7b) 

'combining equations (7b) and (6b) then letting MathCad 
handle the algebra and plotting, we get the following graph. 



If we examine equation (6a) we see that the form of the 
equation for the xy-plane is that of an ellipse displaced 
some distance from the origin ( 8 ) . 

ORIGINAL PAGE IS 

A*x 2 + B*x + C*y 2 = D OF POOR QUALITY 

or 

A ,# (x + b) 2 + C*y 2 = D (8) 

where A t A T ,B t C,D are constants. 

This form of our equation certainly agrees with the plot 
of the shoreline performed on MathCad as one can easily see 
by inspection. Thus we may conclude that the plot agrees with 

Hr\ 




the equation obtained for the fluid shape. 

One may attempt to dispute the validity of the 
equations; however, recall that the problem of the actual fluid 
sloshing was not solved but rather a simplification of the 
problem. The solution though is still useful in gaining 
insight to the degree of fluid sloshing under the simulated 
conditions that have been chosen. 


n 



COMPUTER SIMULATION OF A POINT MASS 


A theoretical analysis representing a -fluid by a 
displaced mass point has been developed by reference 5. 
Using this knowledge to validate our efforts, it was desired 
to track the motion of a point mass in a spinning and 
precessing hemispherical container with a nutation angle 
equal to zero. This was done in reference 6, (see appendix 
S), resulting in the following governing equations: 

- (^) l 5\hdCos6 + 4 (9* 

} t 

= J- \ - IBSMb +k$r - a«_ ) - l) % ] 4 

it % Si re l ” r ' 

l ^ =lLi 

c k Tt 


where 


Using Heun's method -for solving ordinary differential 
equations, as in reference 6, and noting that the quantities 
g , and a are near 1 v constant, the position of the point 
mass as a function of time was *ound. A computer program 
was written to generate the aata needed to develop a plot of 
the motions ( se e appendix C ) . 

Several nitial conditions were used ' , $ . 20, 40, and 

30 RPMs; and friction factors of 0.50, 0.75, 1.50) ana plots 


were m a a e for each of 


*■ K 


e nine con a i t : ons . 


It was f oun< 


that no matter what value the friction factor had, the point 
mass flew out of the cup at SO RPM. In £ act, oarely a 
difference ~ motion was observed (see figures 1-S>. 

However, at 40 RPM, *:he = *;ze of the rr l ctlcn factor 
nfluenced motion in two distinct 'ways, see ~ ' gures 4-6. 
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At a lower value, the motion rose higher in the cup as 
expected. Also as expected because of the greater motion, 
the point mass in the case of lower friction factor took 
much 1 onger to stabi 1 i ze and have a smooth mot i on . 

In the case of spinning at 20 RPM, see figures 7-9, 
much less is obvious. Since the spin rate is so low, the 
point mass stays much closer to the bottom of the cup. The 
most interesting observation is that again the friction 
factor plays a large role in compactness of the motion. 

This theoretical analysis, though it is valid, will not 
have much of an influence in the report except that it shows 
that a stable motion can be achieved. Applying this to our 
problem we can hope that an equilibrium of the fluid 
motion can be found. 
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HOT WIRE 


In the process of designing the experimental fluid 
sloshing machine, the question of how to measure the 
depth of fluid in the sloshing container turned out to 
be a very challenging project. Professor Meyer proposed 
the use of a heated tungsten wire inside of the fluid 
container (cup). The heat loss b$ the wire would in 
turn be related to the depth of fluid inside of the cup. 
Other proposals, like pressure tabs and simple capillary 
tubes, were also considered and thoroughly discused. 

Since the hot wire proposal seemed to be the most 
promising path, several calculations were done to relate 
the heat lost the tungsten wire to the actual depth 
of fluid in the cup. 

The hot-wire was intended to be an application of the 
Hot-Wire Anemometer commonly used in measuring the speed 
of fluids. The hot wire anemometer is basically a 
thermal transducer. An electric current is passed 
through a fine filament which is exposed to a cross 
flow. (This fine wire is actually one of the 
resistances in a Wheatstone bridge circuit). As the 
flow rate varies, the heat transfer from the wire to the 
flowing fluid varies (increases with increasing velocity 
and decreases with decreasing velocity). This variation 
occurs because the electrical resistivity of the wire is 
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a very strong function of temperature. Hence, when the 
wire loses heat (it cools) its electric resistivity goes 
down (true for metals). 

There are basically two techniques to monitor flow 
conditions with the hot wire: constant temperature and 
constant current flowing through the wire. When the 
current in the wire is kept constant, the changes in its 
electrical resistivity unbalance the Wheastone Bridge. 
This is recorded as a voltage drop across the bridge. 
On the other hand, when the temperature of the wire is 
kept constant, then a feed-back control will have to be 
part of the Wheatstone bridge. This feed-back control 
will sense the increase or decrease of heat transferred 
by the wire to its surroundings, and will adjust the 
amount of current flowing through the wire in order to 
keep the temperature constant. 

The problem with the hot wire proposal is not one but 
many. For example, it is necessary to determine whether 
the wire is losing heat through a free convection 
process or a forced convection process. Moreover, it is 
critical to determine how much heat is lost to the 
sloshing fluid, and how much is lost to the surrounding 
air. This is important because if the difference 
between the heat lost to air and to water is not very 
significant then the recorded datum would be misleading. 





4V, eridj* 




Feed Bock, 
Power Sovtce, 


Amplifier. 



Fig. A Block diagram of a constant temperature 
anemometer. The hot wire is the probe 
acting as one of the resistors in the 
Wheatstone Bridge circuit. The feed back 
control adjusts the current to keep the 
bridge balanced. 
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Fig. B The block diagram of a constant current 
anemometer. The probe is one on the 
resistors in the Wheatstone Bridge circuit. 
The voltage across the bridge is shaped and 
amplified before being recorded. 












The following calculations will show that, in view of 
all the "forced" assumptions to idealize the process, 
the hot wire is not the most reliable method to 
determine the depth of fluid in the cup. 
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Now if the same wire were exposed only to air, the air 
could not be assumed stagnant with respect to the wire 
since the cup is rotating about a point ten inches away, 
(see fig.C) 



Hence, the velocity of the air hitting the vertical wire 
is approximately given by the tangential velocity of the 
wire: 


Vqf= R-U>, UA) 

However, this velocity indicated by eg. ( ) does not 

take into account the fact that the wire is also 
rotating about its own axis. So, the absolute velocity 
of a point on the surface of the wire is the sum of V 
given by eq. ( ) and the velocity of the point with 
respect to the 


center 


of 


the 


wire : 
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This ratio is very likely to go down when more of the 
wire is exposed to air as it would be the case during 3° 


actual experiment. So, 
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In view of all the assumptions and of the given sample 
ratios, this proposal of the hot wire presents too many 
uncertainties. For this reason no experimental work was 
carried on and the proposal had to be turned down. 
However, there is another proposal that is likely to 
give more precise results in a simpler manner. This 
proposal, given by James Marcolesco (TA), will be 
explored further this Spring. 






The? purpose of these calculations is to demonstrate that 


the heat loss of tungsten wire is greater in water than in 
air and hence the voltage drop is also greater in water. 
However * a very simplified version of calculating the heat 
transfer of the hot-wire probe is as follows: 

A tungsten wire of length * L » 7.5mm and diameter? D - 1.27x 
lO^mm is used. There are two cases where this tungsten wire 
is either exposed to the air or water in the rotating cup 
with no nutation angle. Several assumptions are made due to 
the complexity of these problems? such as forced convection 
heat transfer of wire? the wire is entirely exposed to the 
air or submerged in water ? the temperature of wire is below 
the boiling-point of water? and both wire and cup have equal 
angular velocity (or constant velocity) on the stationary 
arm . 


The sample calculations of these two cases are as 
fol lows: 

I. To find the velocity of water in the cup using Navier- 
Stofces equation. 
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We can observe that the heat loss of wire in the air is 
much less compared to the heat loss of wire in the water. It 
can therefore be argued that the heat loss in air is Quite 
nea 1 i a i b 1 e . 

After performing these sample calculations? we still 
need to examine the effects of heat loss, and also the voltaae 
difference of tungsten wire having a different height 
submerged in water. Let us make several assumptions 1 oi this 
case bv: 


1. neglecting heat loss of tungsten wire in the air (as 

shown earlier where heat loss in air is much smaller than 
in water ) . 

E. havina the temperature of wire below the boiling-point of 
water . 

3. having constant velocity (at V 0 = 0.63m/s). 

We have the same conditions for wire and water 
temperatures? except that the height of wire submerged in 
water is 5.0mm. This sample calculation follows: 


. Forced convection of wire in water, 

5-1 


h\ r 
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Exaraininn these three cases* we can arque that for 
smaller values of tungsten wire height being submerged in 
water, we will obtain lower heat transfer coefficients, lower 
values of heat loss and thus, lower voltage differences. So, 
based on these "rough" assumptions that as the convection 
heat transfer coefficient decreases, so will the voltage 
differences, and we can try to correlate these to give h = 
h (V) . 
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CONCLUSION 


Although three-axis stabilization presents -few sloshing 
effects, the ease of spin stabilization means that liquid 
sloshing still remains a problem. Nutational stability 
characteristics have generally been established by costly, 
■full scale tests. By use of the "follower force" method, 
results will be more easily obtained and just as accurate. 
Two different simplified methods for the observation of 
liquid slosh will be used in this design. Pictures will be 
taken fc- qualitative analysis and the depth of fluid at 
various points in the cup will be used for quantitative 
anal ysis. It is hoped that after the mode 1 i s bu i 1 t , the 
results will be helpful in gaining an insight into the 
problem of liquid sloshing. 
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Parameters of the modeled system 


O := 0,-08 . . 2 • tt ra := 15 cm 
CONST := .91 cm 0 := 0 


rc := 15 -cm Me := 80‘rpm wa := ^Orpm 
R : = rc 


"The equation representing the fluids' shape: 


L ( R , 0 , 0 ) 



2 2 

Wa + Me 

ra + R • cos (0) • cos ( 0 ) 

2 


■ R • c o s ( 0 ) • c o s ( 0 ) 


2 

2 Me 

+ (R • cos (0 ) ' sin(0) ) + CONST 

2 ' g 


"The equation of a spnere: 
P ( R , 0 , 0 ) : = R 


"Spherer ical Coordinates : " 

: : ( R , 0 , 0 ) = R cos(0) cos(O) 
v ( R , 0 , 0 ) = R • cos ( 0 ) • s in < 0 ) 
z ( R , 0 , 0 ) = R s i n ( 0 ) 




<,?, o ) , :■ 


s ,e>,x\P<R,*,0>,Pf®> 15-cm 
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program bal 1 

C*************** MAIN PROGRAM 
C*************** 

implicit double precision (a-z) 

parameter (pi = 3. 141592654, neqns=4,eps = 1 .00-06) 

common wr , w the t , wps i , w the 1 1 , wps it,gthet,gpsi ,athet,aps 
dimension LK10) 

integer n , nsteps ,max i t , i , j , np i 


open< 15, f i 1 e = 

'xy. tgp') 

open ( 20 , f i 1 e = 

"xyz .out' ) 

open < 1 6 , f i 

1 e = 

* input. dat y ) 

read<16,*> 

rpm 


r ead< 16,*) 

r 1 


r ead< 16,*) 

r 


read< 16,*) 

r0 


read< 16,*) 

r 2 


read< 16,*) 

r 3 


read< 16,*) 

nu 


read( 16,*) 

u(l) 


read( 16,*) 

u ( 2) 


read< 16,*) 

u < 3) 


read< 1 6 , *) 

u ( 4) 


read< 1 6 , *) 

ta 


read<16,*) 

tb 


print * , " 

input 

max i t and h 

read * , max i t , h 



nsteps = abs< ( tb-ta)/h ) 
wl=rpm*2 . 0*p i /60 . 0 

do 10, n=Q .nsteps 
t=tb*n/nsteps 

print*, ' n , u < 1 ) , u ( 3) : 7 , n , u < 1 ) , u < 3) 

call Params( wl , r , r 0 , r 1 , r2 , r 3 , U , t ) 
call Heun < neqns ,max i t , h , eps , t , U) 

the = LKl) 
psi = U( 3 ) 

if ( dabs ( t he ) . gt . 2 . 0*p i ) then 

np i = di n t< the/2 . 3/p i ) * 

the = the-np i *2 . 0*p i 
end i f 

if < dabs< p s i ) . g t . 2 . 0*p i ) then 
npi = d i n t ( ps i /2 . 0/p i ) 
psi = ps i -np i *2 . 0*p i 
end i f 

LKl) = the 
IK 3) = psi 


x = r*sin<U<l))*cos<U<3)) 
y = r*si n<U< 1 ) )*si n<U<3) ) 
z = r*cos( IK 1 ) ) 
if <MOD<n ,50) .eq .0) then 
write <15,1000) x,y 
write <20,1002) n,t,x,y,z 
1000 format < e 1 5 . 8 , 2x , e 1 5 . 3) 

1002 format < i 6 , 2x , e 1 0 . 5 , 2x , e 1 5 . 8, 2x ,el 5.8 , 2x , e 1 5. 8) 

endi f 

if <U< 1 ) .gt . <p i/2.0) ) then 
stop 
endi f 

if <U< 1 ). 1 t . <-p i/2.0) ) then 
stop 
end i f 
10 continue 
end 


C*********** THE VECTOR FUNCTION 
C*********** 

double precision function f<i,U,t) 
implicit double precision < a-z ) 

common wr,wthet,wpsi ,wthett,wpsit,gthet,gpsi ,athet,apsi 
i nteger i , np i 
dimension U<10) 

print *, ' f unc ' ,U<1) ,U<3) 

if < i . eq . 1 ) then 
f = U< 2) 

else if (i.eq.2) then 

f = U< 4) *U< 4) *s i n < U< 1 ) ) *cos< U< 1 ) ) + <gthet - athet)/r 
% nu*< 2) - wr*<wthe t - 2 . 0*U< 4) *s i n< U< 1 ) ) ) - wthett 

else if <i.eq.3) then 
f = U< 4) 

else if < i . eq .4) then 

f = <-2.0*U<2)*U<4)*cos<U< 1 ) ) + <gpsi - apsi)/r - 
y. nu*U<4)*si n<U< 1 > ) -wr*<wpsi + 2.0*U<2)> 

y. + wthett)/si n<U< 1 ) ) 

end i f 
end 


C********** THE FORMULAS FOR THE PARAMETERS 
C********** 

subr ou tine Par ams< wl ,r,r0,rl , r-2 , r 3 , U , t ) 
implicit double precision <a-z) 

common wr , w the t , wp s i , wthett, wp si t,gthet,gpsi , a the t , aps i 
integer np i 
dimension U<10) 


print *, ' par am' , U( 1 ) , U< 3 ) 


9 

hat = 
mx = 
wy = 
mz = 
wr = 
X 

wthe t 

x 

wps i 

MX t 
My t 
MZ t 

Mthe 1 1 

X 

wpsi t 
9 * 

9 / 

9 Z 

gthe t 

gpsi 

ax 

ay 

az 

athe t 

X 

aps i 
end 


9.81 

a t an < m 1 *m 1 *r 1 /g ) 

Ml*sin<U<l) ) *cos<r2*Ml *t/r3) 

Ml*sin<U<l))*sin < r 2 *m 1 *t/r3) 

-Ml*<r2/r3+cos<U<l>>> 

MX*si n<U< 1 > )*cos<U< 3 ) ) + My*s i n < U< 1 ) ) *s i n < U< 3 ) ) + 
MZ*COS< U< 1 ) > 

= MX*cos(U< 1 ) >*cos<U<3> ) + My*cos( U< l))*sin<U<3)) - 
wz#s i n ( U( 1 ) ) 

= -wx*s i n < U< 3> ) + My*cos<U<3>> 

= -Ml*Ml*r 2 /r 3 *si n<hat>*si n ( r 2 *Mi *t/r 3 > 

* m1*m 1 *r2/r3*si n < hat ) *cos< t2*m 1 *t/r3> 

= 0.0 

= Mxt#cos<U< 1 ) >#cos<U<3> ) + Myt*cos<U< 1 ) > *s i n<U<3> > - 
mz t*s i n < U< 1 > ) 

= -mx t*s i n < U< 3) ) + My t*cos( U( 3) ) 

= -g*s i n ( hat ) *cos< t2*m 1 »t/r3) 

= -g«s i n < hat ) *s i n < t2*m 1 *t/r3) 

=s g*cos<hat) 

= gx*cos< U< 1 > ) *cos< U<3) > + gy*cos<U< 1 ) >*si n<U<3) ) - 

gz*s i n < U( 1 ) ) 

= -gx*s i n ( U( 3> > + gy*cos<U<3)> 

ss -Ml *m 1 *r0*cos( U< 1 > >*cos<r2*Ml*t/r3) 

= -Ml * m 1 *r0#cos( U< 1 ) ) *s i n < t2*m 1 *t/r3) 

= -Mi*Ml*rO*sin(U<l>> 

= ax»cos(U< 1 ) >*cos<U<3) ) + ay*cos< U< 1 ) ) *s i n< U< 3) ) 
-az*s i n < U< 1 > ) 

= -ax*s i n < LK 3) ) + ay*cos<U<3>> 


C********** SUBROUTINE JACOBN 
C***«****»* 


subroutine Jacobn < J , U, t ) 

common wr ,Mthe t ,wps i ,Mthett,Mpsit,gthet,gpsi , athe t , aps i 
dimension J(10,10) 
d i mens ion U( 1 0 ) 

implicit double precision (a-z) 
integer npi 

print *, ■'jaco', U(1),U(3) 


J ( 1 , 1 > 
J< 1 ,2) 
J<1 ,3) 
J< 1 ,4) 
J( 2 , 1 5 
J<2, 1 ) 
J< 2 , 2) 
J<2,3> 
J ( 2 f 4) 
J < 3 , 1 > 
J ( 3 , 2 ) 


0 

0 

0 

0 

U<4)*U(4>*<cos<U< 1 ) >*cos< U< 1 ) )-s i n(UU ) >*si n<U( 1 > ) ) 


J<2, l)+wr#( 2*U< 4)*cos<u<l ) ) ) 
-NU 


0 

2#U< 4) «SIM< U< 1 ) > *C0S< U< 1 > >+2*Mr*s i n<U< 1 ) ) 


0 

0 




J< 3 , 3) = 0 
J< 3 , 4) = 1 

J< 4 , 1 ) - -2*U< 2) *U(4>*cos<LKl)) + <<gpsi -aps i )/r ) 

J< 4 , 1 ) = J<4,l)-nu*U(4)*sin<U(l)> -wr*(wps i + 2*U< 2))+wthett 
J<4, 1 ) - J< 4 , 1 ) *cos( LK 1 ) ) 

j< 4 , i j n (U< 1 ) )*<2*U<2)*U<4>*si n<LK 1 ) ) -nu*U< 4) *cos< U< 1 > > )-J<4, 1 ) 
J( 4 , 1 ) = JC4 v l>/<sin(U(l))*«in<U(l))> 

J< 4 , 2) = -2*U<4)*C0S<U< 1 > >/SIN<U< 1 > >-<2*WR/si n<U< 1 > ) > 

J< 4 , 3) = 0 

J<4,4> = < -2*U< 2) *COS< U< 1 ) ) /SINK IK 1 ) ) ) -NU 

RETURN 

END 


n 


c 

c 

c MARBLE : finds the position of a marble placed in the TEACUP 

c 

c James Marcolesco 

c Winter 1990 


c 

MANE 

199 

Pro-f. Meyer 

c 

c 

c 

This 

routine drives the Heun subroutine 

c — 

c 

c 

max i t 


maximum number o-f iterations 

c 

h 

= 

step size 

c 

eps 


convergence tolerance 

c 

neqns 

= 

number of equations to be solved in system 

c 

F 

= 

■function vector 

c 

J 

= 

Jacobian matrix o-f -function vector 

c 

du 

S 

vector to add to solution vector 

c 

U 

= 

solution vector to be solved in HEUN 

c 

c 

c 

t 


current time at call to HEUN 


c 234 wi- 7 


progr am ma i n 

implicit double precision (a-z) 

integer max it, neqns, nsteps,i , j , n ,np i ,r ad, count,! imit 
common r , nu , gthe , gph i ,athe,aphi , w3r , wthe , wph i , wph i dt ,wthedt 
dimension u < 1 0 ) 


c 


c 

c 

c 


open <15, -file = "xy.tgp" 

open <16, -file = "x.out" 

open <17, -file = '/.out' 

open <18, -file = "z.out" 
open <19, -file = "xyz.out" 

open <20, -file = "input.dat 


, status = "unknown' ) 

, status = "unknown") 

, status — "unknown") 

, status = "unknown") 

, status = "unknown") 

, status = "old") 

. . . . .set control parameters 


neqns = 4 
g = 9 . 80 665d0 

pi = 3.141 592654d0 

p i 2 = 2 . dO * p i 

p j 90 = pi / 2 . dO 
rtd = ISO.dO / pi 


wr 

pr 

pr 

wr 

pr 

pr 

wr 


te<*, 1002) 

n t* , "Enter max iterations:" 


te<*, 1001 ) 

r<t*, "max I t *» 1 is explicit Heun method" 

nt*, "max it > 1 is implicit trapezoidal method" 

t e < * , 1 0 0 1 ) 


r e ao+ , max ; * 





wri te(*,1002) 

print*, 'input convergence tolerance, eps:' 
wri te(*, 1001) 
read*, eps 

write (19,1002) 

i -f (max i t . eq . 1 ) write<19,*) '** Explicit Heun-s Method **' 

1 4 (mtxlt.gt.l) write<l? f *> ' ** Implicit Trapezoidal Method **' 
write (19,1002) 

1001 format (/) 

1002 format (//) 

1 c on t i n u e 
count = 0 
wri te ( * , 1002) 

print*, 'input time step, h:' 

wri te ( * , 1001) 

read*, h 

wri te<*, 1002) 

print*, 'for plot files, skip how many points?' 
wri te ( * , 1001) 
read* , limit 
wri te(*,1002) 

read( 20 , *) RPM 
read( 20 , *) r 1 

wl = RPM * 0 . 1 0472d0 

THETA = datan ( wl * wl * rl / g) 


read( 20 , *) r 
read( 20 , *> R0 
read( 20 , *) R2 
read( 20 , * > R3 


c 

c 

c 


re ad (20,*) nu 


read in initial conditions 


do 2 , i = 1,4 
read(20,*) u(i) 
continue 


read(20 , *) ta 
read(20,*) tb 

the = u ( 1 ) 
phi = u < 3 ) 


wr 


19,1000) h ,THETA*rtd,wl , r , r 1 , R0 , R2 , R3 , nu , t h e , p h i 


10 00 format : lx, step size, h = •' ,f 7.4, sec 

% 
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, f 7 . 4 , 7 rad/s 


c 

c 

c 


c 

c 

c 

c 

c 


1 x , 'THETA = ',i7.4,' deg^Sx^wl » ' 

2 x/r =',**• 4,' m^/Zrl -',*6.4,' m' , 

. 2x f / RO =',f6.4,' m / , 1 x , / R2 -',*6.4,' m' 

2x, / R3 =',*6.4,' m',2x,/,'nu =',*6.4, 
2x,'the<0> =' ,f6.4, ' rad" , 2x , ' ph i ( 0 ) 

. 3x , ' t ime step ' ,5x 9 ' t i me ' , ?x , ' x ' 9 12x , ' y' 




,*6.4,' rad' 
, 14x , 'z' ) 




! 


(1) calc nsteps 


nsteps = dabs(tb-ta) / h + 1 
t = ta - h 


(2) begin time stepping 
include t = ta 
Calculate coe+*icients... 


do 100, n ** 1 , nsteps 

t = t + h 

ms in = ds i n < the ) 
mcos = dcos<the> 

ps i n = ds i n< ph i ) 
pcos = dcos( ph i ) 

tsin = dsin<THETA> 
tcos = dcos<THETA) 

rsin = dsin<R2 / R3 * wl * t) 
rcos = dcost R2 / R3 * wl * t) 

w3x = wl * tsin * rcos 
w3y = wl * tsin * rsin 
w3z — — wl * ( R2 / R3 + tcos) 


w3r 

= Vvi3x * ms i n 

* 

pcos + 

w3y * ms i n 

* psin + w3z 

w3the 

= w3x * mcos 

* 

pcos + 

w3y * mcos 

* psin - w3z 

w3ph i 

= - w3x * psi 

; n 

+ w3y 

* pcos 



w3xdt = - wl * wl * R2 / R3 * tsin * rsin 
w3ydt * wl * wl * R2 / R3 * tsin * rcos 
w3zdt = 0 . dO 

wthedt = w3xdt*mcos*pcos + w3ydt*mcos*ps i n - w3zdt*msin 
wphidt = - w3xdt * psin + w3ydt * pcos 


gx = - g * tsin * rcos 
gy = - g * tsin * rsin 
gz = g * tcos 


gthe “ gx * mcos * pcos + gy * mcos * psin - gz * msin 
gphi = - gx * psin + gy * pcos 


a ux = - wi 



Ai i 
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aOy = - w 1 * wl * rO * tcos * rsin 

aOz = - wl * wl * rO * tsin 

aOthe = aOx * mcos * pcos + aOy * mcos * psin - aOz * msin 

aOphi = - aOx * psin + aOy * pcos 

c 

call solving subroutine 

c 

call heun<neqns,max i t ,h ,eps, t ,U) 

the = u<l) 
phi = u < 3) 
c 

c . reduce to 3rderk2 pi) 

c 

if ( dabs< the ) . gt . p i 2) then 
np i = d i n t < the/p i 2) 
the = the - npi * pi2 
end i f 


if ( dabs< ph i ) . gt . p i 2) then 
npi = dint<phi/pi2) 
phi = phi -npi *pi2 
end i f . 

find x,y,z coords 

x = r * dsin<the) * dcos<phi) 
y = r * dsin<the) * dsin(phi) 
z = r * dcos( the ) 

.....(11) print out results 


if (count.gt. limit) then 
wr i t e ( 1 5 , 1 1 0 2 ) x,y 

write <16, 1101) x 
write (17, 1101) y 
wr ite<18,1101) z 

write<l?,1105> n,t,x,y,z 
count = 0 
end i f 

1101 f ormat < 1 x , e 1 5 . 3) 

1102 f ormat < 1 x , e 1 5 . 3 , 2x , e 1 5 . 8) 

1105 format < 4x , i 5 , 7x , f 6 . 3 , 3< 2x , e 1 2 . 4) ) 
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stop if out of teacup 


if < dabs( the ) . gt . p i 90 ) then 
wr i te < * , 1 000 ) 

print*, mar cl e flew cut of teacup: "META = 


, th e*r td , " deg-' 




c 

c 

c 


stop 
end i f 

u(l) = the 
u ( 3) = phi 

100 continue 

c 1 ose (15) 
c 1 ose < 1 6) 
c 1 ose (17) 
c 1 ose (18) 
c 1 ose < 1 ?) 
c 1 ose < 20 ) 

end 


c functions 


double precision function f(i,U,t) 
implicit real*8(a-z) 
integer i 

common r , nu , gthe , gph i ,athe,aphi , w3r , wthe , wph i , wph i dt , wthedt 
dimension u ( 1 0 ) 

cl = dcos( u(i>) 
c2 = cl * cl 
si = dsin(u(l>) 
s2 = si * si 

goto <10,20,30,40), i 

10 f » u<2) 
return 

20 f = u<4> * u<4) * cl * si + (gthe - athe) / r - nu * u<2) 

- w3r * (wthe - 2.d0 * u<4) * si) - wthedt 
return 

30 f = u(4) 
re turn 

40 f = ( -2 . dO * u<2) * u(4) * cl + < gph i -aph i )/r - nu * u<4) * si 
. — w3r * ( wph i + 2.d0 * u(2)) + wthedt) / si 

return 

end ORIGfNAL PAGE fS 

OF POOR QUALITY 

c 

c Jacobian subroutine 


subroutine Jacobn ( J , 'J , t ) 
implicit real*8(a-z> 

; ommon r , n u , g t h e , gp h i , a t h e , ap n : , won , w t n e , wp n ! < wp n i ct , w t h e o •. 


dimension J(10,10>,u<lQ) 


cl = dcos(u<l>> 




c2 = cl 

* cl 





si = ds i. n < u < 1 ) ) 




in 

II 

w 

in 


* si 





j<i ,d 

s 

O.dO 





J<1 ,2) 

= 

1 .dO 





J<1 ,3) 

= 

O.dO 





J< 1 ,4) 

= 

O.dO 





J<2, 1 ) 

= 

u ( 4) 

* 

u < 4) * 

<c2 - s2) + 

2.d0 * u<4) * w3r * cl 

J<2 , 2) 

= 

- nu 





J<2,3) 

= 

O.dO 





J<2 , 4) 

= 

2 . dO 

¥r 

u < 4) * 

cl * si + 2. 

dO * w3r * si 

J< 3 , 1 ) 

= 

O.dO 





J< 3 , 2) 

= 

O.dO 





J< 3 , 3) 


O.dO 





J<3,4) 

= 

1 .dO 





J<4, 1 > 

= 

<2 

.dO 

* u ( 2) 

* u < 4) * si 

- nu * u (4) * c 1 > / si 



- c l/s2 

* (-2. 

dO * u<2> * 

u < 4) * cl + (gphi - aph i ) / r 





— nu * u ( 4) * si — w3r * (wphi + 2.d0 * u(2)) 





+ wthedt) 


J< 4 , 2) 

= 

< -2 . 

dO 

* u ( 4) 

* cl - 2 . dO 

* w3r> / si 

J< 4 , 3) 

= 

O.dO 





J( 4 , 4) 

=r 

<-2. 

dO 

* u < 2) 

* cl - nu 

* si ) / si 


return 

end 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 


HUEN's METHOD 


James Marcolesco 
Fall 1989 

MANE 1 92C Prof. McDonough 

This routine solves Initial Ual ue Problems -for ODE'S 


m 

max i t 
h 

a , b 


number of final iterations 
maximum number of iterations 
step size 
t i me doma i n 


F 

J 

du 

U 


function vector 

Jacobian matrix of function vector 
vector to add to solution vector 
solution vector to be solved in NEWTON 


±m. W I 


V 
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subroutine heun(neqns .max i t,h,eps,t,U> 
implicit double precision <a-h,j,o-z) 


dimension Jf<10,10), JFF <10,10), FF< 10) ,u(10) , uol d< 1 0 ) , du < 1 0 ) , 

. ustar<10) , g< 10 ) 

double precision maxdif 

integer de 1 ta , i , j ,m ,max it,n,neqns,nsteps 
external f ,de 1 ta 
c 

.<3) begin Newton iterations 

c 

3 tol d = t - h 
m = 0 

do 80 , m = 1 , max i t 
if <m . gt . 1 ) goto 6 
c 

< 4) Evaluate Ustar -for use in 

c Huen's method 

4 do 20 , ' = 1 , neqns 

g(i) = f ( i , u , tol d) 
ustar(i) * u<i) + h * g< i ) 

20 continue 

c 

C <5) Calculate inital guess -For 

c trapezoidal rule -from Huen's 

c method 

c 

5 do 30, i = 1 , neqns 

uold(i) = u<i) 

u<i) = u(i> + 0.5 * h * (g<i) + f<i,ustar,t)) 

30 continue 

i-f Cmaxit.eq.l) return 
c 

< 6 > Load J < f ') -or ' Jew * on t 5 r a : 

c 

6 call Jacobn < J*f , U , t ) 
c 

.............. <7> Evaluate FF ( Um ) and J ( F ) 

c 

7 do 50, i = 1 , neqns 

FF< i ) = u ( i ) - 0.5*h*f < i ,u,t) - <uold<i) + 0.5*h*g(i)> 
do 40 , j = 1 , neqns 

,JFF<i,j) = de 1 ta< i , j ) - 0.5 * h * Jf<i,j> 

40 continue 

50 continue 

c 

<3) Solve -for du using Gauss 

c 

3 call GAUSS (neqns , JFF , du , FF) 

..(9) Calculate max norm of du< i ) 

an c : n cr erne nr 


ORIGINAL PAGE IS 
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9 dumax = 0 . dO 

do 60 , i = 1 , neqns 

i f ( dabs( du < i )>. gt . dumax ) dumax = dabs<du<i>) 
u< i > = u<i> - du< i > 

60 continue 


c 

< 10) Test convergence of Newton 

c iterations 

c 

10 if (dumax.lt.eps) return 

80 continue 


print*, 'Newton iterations failed to converge at time step n + 1' 

return 

end 

c 

c define kroniker delta function 

c 

f unc t i on del ta< i , j ) 
implicit integer(a-z) 


if < i . eq . j ) delta = 1 
if < i ,ne . j ) delta = 0 


return 

end 


c GAUSS ELIMINATION SUBROUTINE 

c 

c James Marcolesco 

c MANE 1 92C 

c Fall 198? Prof. McDonough 

c 

c 

C234567 

subroutine gauss(n,A,X,B> 
implicit rea 1*8 a-h , m , o-z ) 
dimension A<10,10),X<10),B<10),M<10,10) 
c 

c ......... . forward el i m i nat i on 

c 

do 100, k = 1, n-1 


c ............ r ow pivoting 

r 

imax = K 

amax = dabs( A< k , k > > 


% 


do 10 , i = k + 1 , n 

i i < dabs< A< i , k ) ) . gt . amax ) then 
amax = dabs<A<i,k>> 
i max = i 
end i f 

10 continue 

i i (imax.eq.k) go to 30 
do 20 , j * k , n 
a temp = A< k , j ) 

A<k , j ) = A< imax , j ) 

A< i max , j ) = atemp 
20 continue 

btemp = b<k) 
b< k ) = b( i max ) 
b< i max ) = btemp 

30 do 60 , i = k + 1 , n 

M< i ,k) = A< i ,k> / A<k ,k) 

b< i ) » b< i ) - M< i ,k) * b<k> 

do 40 , j = k + 1 , n 

A< i., j) = a< i ,j> - M< i ,k) * A<k,j> 
40 continue 

60 continue 

100 continue 

back substitution (solve stage) 

x(n) = b(n) / A(n,n) 

do 200, i = n - 1, 1, -1 

x( i ) = 0 .d0 

do 21 0 , j = i + 1 , n 

;<(i) = x < i ) + A < i , j ) * x(j) 

210 continue 

x(i) * (b(i> - x< i >) / A<i,i) 

200 continue 

return 

end 
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APPENDIX e 


Calculation of velocity of sloshing fluid around wire 

Assume two concentric cylinders. The inner cylinder is 
fixed to the outer one, and both of them are rotating 
with the same angular velocity about their common axis 
of symmetry. Moreover, assume that the velocity of the 
fluid around the inner cylinder is only a function of 
the distance r from the center of this cylinder 
(wire). For simplicity no flow in the direction of the 
axis of symmetry is assumed (this is noi* a very good 
assumption), and the cylinders are very long in their 
axial direction. Also, as it will be shown, the main 
problem with this analysis is the Boundary Conditions. 
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